
function [put,delta] = BS_put(S0,K,r,T,sig)

d1    = ( log(S0/K)+(r+sig^2/2)*T ) / ( sig*sqrt(T) );
d2    = d1 - sig*sqrt(T);
delta = normcdf(d1)-1; % = -N(-d1)
put   = K*exp(-r*T)*normcdf(-d2) + delta*S0;

end
